﻿
using System;

namespace EmperialApps.WeatherSpark.Data {

    /// <summary>Calculates astronomical positions and phases.</summary>
    public static class Ephemerides {

        /// <summary>Gets the Julian date for the specified local date.</summary>
        public static double GetJulianDate( DateTimeOffset date ) {
            DateTimeOffset offsetDate = date.AddDays( GregorianJulianOffset )
                .ToUniversalTime( );
            TimeSpan startDifference = offsetDate - JulianStart;
            return startDifference.TotalDays;
        }

        /// <summary>Gets the time of sunrise and sunset for the specified date.</summary>
        /// <remarks>Implementation of the sunrise equation.</remarks>
        /// <see href="http://en.wikipedia.org/wiki/Sunrise_equation"/>
        /// <see href="http://www.esrl.noaa.gov/gmd/grad/solcalc/sunrise.html"/>
        public static Pair<DateTimeOffset, DateTimeOffset> CalculateDaylight( Coordinate location, DateTimeOffset date ) {
            double latitude = location.Latitude;
            double longitude = location.Longitude;
            double timezoneOffset = date.Offset.TotalHours;

            double julianDay = GetJulianDate( date ) + 0.5;
            double julianCentury = (julianDay - 2451545) / 36525;

            double solarMeanLongitude = (280.46646 + julianCentury * (36000.76983 + julianCentury * 0.0003032)) % ConvertValue.DegreesPerCircle;
            double solarMeanAnomaly = 357.52911 + julianCentury * (35999.05029 - 0.0001537 * julianCentury);
            double earthOrbitEccentricity = 0.016708634 - julianCentury * (0.000042037 + 0.0000001267 * julianCentury);
            double solarEquationOfCenter =
                      Sine( solarMeanAnomaly ) * (1.914602 - julianCentury * (0.004817 + 0.000014 * julianCentury))
                + Sine( 2 * solarMeanAnomaly ) * (0.019993 - 0.000101 * julianCentury)
                + Sine( 3 * solarMeanAnomaly ) * 0.000289;

            double solarTrueLongitude = solarMeanLongitude + solarEquationOfCenter;
            double solarApparentLongitude = solarTrueLongitude - 0.00569 - 0.00478 * Sine( 125.04 - 1934.136 * julianCentury );
            double meanObliqueEcliptic = 23 + (26 + ((21.448 - julianCentury * (46.815 + julianCentury * (0.00059 - julianCentury * 0.001813)))) / 60) / 60;
            double obliqueCorrection = meanObliqueEcliptic + 0.00256 * Cosine( 125.04 - 1934.136 * julianCentury );

            double solarDeclination = ArcSine( Sine( obliqueCorrection ) * Sine( solarApparentLongitude ) );
            double y = Math.Pow( Tangent( obliqueCorrection / 2 ), 2 );
            double equationOfTime =
                  4 * ConvertValue.ToDegrees( y * Sine( 2 * solarMeanLongitude )
                - 2 * earthOrbitEccentricity * Sine( solarMeanAnomaly )
                + 4 * earthOrbitEccentricity * Sine( solarMeanAnomaly ) * y * Cosine( 2 * solarMeanLongitude )
                - 0.5 * y * y * Sine( 4 * solarMeanLongitude )
                - 1.25 * earthOrbitEccentricity * earthOrbitEccentricity * Sine( 2 * solarMeanAnomaly ) );
            double sunriseHourAngle = ArcCosine( Cosine( 90.833 ) / (Cosine( latitude ) * Cosine( solarDeclination )) - Tangent( latitude ) * Tangent( solarDeclination ) );

            double solarNoon = (720 - 4 * longitude - equationOfTime + timezoneOffset * 60) / MinutesPerDay;
            double sunrise = (solarNoon * MinutesPerDay - sunriseHourAngle * 4) / MinutesPerDay;
            double sunset = (solarNoon * MinutesPerDay + sunriseHourAngle * 4) / MinutesPerDay;

            var daylight = Pair.Create(
                date.AddDays( sunrise ),
                date.AddDays( sunset ) );
            return daylight;
        }


        #region Private Members

        private const double MinutesPerDay = 1440;
        private const double GregorianJulianOffset = 2415018.5;
        private static readonly DateTimeOffset JulianStart = new DateTimeOffset( ).AddDays( 693593 );


        private static double Sine( double degrees ) {
            double radians = ConvertValue.ToRadians( degrees );
            return Math.Sin( radians );
        }

        private static double Cosine( double degrees ) {
            double radians = ConvertValue.ToRadians( degrees );
            return Math.Cos( radians );
        }

        private static double Tangent( double degrees ) {
            double radians = ConvertValue.ToRadians( degrees );
            return Math.Tan( radians );
        }

        private static double ArcSine( double sine ) {
            double radians = Math.Asin( sine );
            return ConvertValue.ToDegrees( radians );
        }

        private static double ArcCosine( double cosine ) {
            double constrained =
                cosine > +1.0 ? +1.0 :
                cosine < -1.0 ? -1.0 :
                cosine;
            double radians = Math.Acos( constrained );
            return ConvertValue.ToDegrees( radians );
        }

        #endregion
    }

}
